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Abstract 

The 0^ real scalar field theory on a fuzzy sphere is studied numerically. We refine the phase 
diagram for this model where three distinct phases are known to exist: a uniformly ordered phase, a 
disordered phase, and a non-uniform ordered phase where the spatial 50(3) symmetry of the round 
sphere is spontaneously broken and which has no classical equivalent. The three coexistence lines 
between these phases, which meet at a triple point, are carefully located with particular attention paid 
to the one between the two ordered phases and the triple point itself. In the neighbourhood of the triple 
point all phase boundaries are well approximated by straight lines which, surprisingly, have the same 
scaling. We argue that unless an additional term is added to enhance the effect of the kinetic term the 
infinite matrix limit of this model will not correspond to a real scalar field on the commutative sphere 
or plane. 
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1 Introduction 



The fuzzy approximation scheme lUl consists in approximating the algebra of functions on a manifold 
with a finite dimensional matrix algebra and in principle provides a regularization of field theory on 
this space which can be used as an alternative to discretising the underlying space as is done in lattice 
field theory. Both the two-dimensional commutative f2^ and Moyal planes can be viewed as the limits 
of a fuzzy sphere of infinite radius. 

Here, we study a real scalar field, (/), with (p^ interaction, in the fuzzy approach using Monte Carlo 
simulations. This becomes a Hermitian matrix model on the fuzzy sphere. The study reveals that 
the model has three distinct phases: (i) A disordered phase; (ii) a uniformly ordered phase and (iii) 
a non-uniformly ordered phase assimilated to a striped phase EHHHl. We find the collapsed phase 
diagram and in particular we calculate the uniform ordered/non-uniform ordered line that was absent 
in |[6l and locate the triple point where the three phases meet. As the mass parameter varies, the non- 
uniformly ordered phase is absent for sufficiently small coupling, but as the coupling is increased this 
new phase opens up between the disordered and uniformly ordered phases. The three phases meet at 
a triple point. 

The transition from the disordered to the non-uniformly ordered phase can also be identified with 
the one-cut to two-cut transition in matrix model theory QIHIll. This transition line merges with 
the predicted curve obtained from the quartic potential of the single trace pure matrix transition for 
sufficiently large couplings, i.e. sufficiently above the triple point. The qualitative features of the 
phase diagram are governed by this triple point. The presence of the non-uniformly ordered phase is 
the principal feature that distinguishes the phase diagram of the fuzzy model from its commutative 
counterpart. 

A preliminary version of these results were presented in Lattice 05 and appeared in lITOl . The 
principal aspects of these results have been confirmed in subsequent studies by Panero iHTl [T2]| and 
Dasetal. l[T3l . 

The current study could be relatively easily repeated for a Hermitian scalar field on other fuzzy 
spaces. The simplest extension would be to fuzzy x or to fuzzy CP^ lfT4l . Fuzzy versions of 

and are also accessible ifTSl and will hopefully be studied in the near future. In all cases, the 
structure of the phase diagram should be similar, although there is no guarantee that all coexistence 
lines will collapse with a consistent scaling as happens for the two dimensional sphere. One prediction 
for the general case is that the disordered non-uniformly ordered line will always be present for 
sufficiently large coupling and will again merge with the pure one-cut two-cut transition for the pure 
matrix model. 
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In section 2 we review the construction of the fuzzy sphere and in section 3 we present the model, 
section 4 describes the Metropolis algorithm, section 5 studies Umiting models such as the lowest ma- 
trix size (two by two) and the pure matrix model, section 6 describes the observables and simulations, 
particularly the specific heat which we use to locate transitions. Section 7 gives our main results and 
describes the collapsed phase diagram and locates the triple point. Section 8 gives our conclusions. 
The paper ends with some technical appendix for the optimization of the simulations. 



2 The fuzzy sphere 

Before introducing the fuzzy sphere, let us look at some basic properties of the ordinary continuum 
2— sphere. A 2-sphere centered on the origin, with radius R, embedded in M^, denoted simply can 
be defined as the set of points (Xj , , Xj ) in M? such that x^^+x^+x^ = R^. It can also be expressed 
by the angles (t?,<p) of spherical coordinates. 

Taking two elements of the algebra, / (jp) and g{'&,q)),vje define their inner product as 

{f\g)= I dar{^,(p)g{^,(p), (1) 

and their norm as 

ll/f = (/!/)= /dn 1/(^,9) 1^ (2) 

where / dH = / d(p / d sin (t?). The norm must be finite for any element of the algebra (square 
integrable functions). Both equations, ([T|l and define the Hilbert space ^ which allows us to 
quantize the theory. 

In general, the Laplace operator contains information on the geometry of the space, i.e. it de- 
pends on the metric as V^- = -j=dj^J\g\d^ where g is the determinant of the metric tensor guv on 
Riemannian and pseudo-Riemannian manifolds |16]. In particular, the Laplacian on the sphere is 

= } + ^2 sin ^ {^^^^7^^- "^^^ eigenfunctions of this operator are the spherical har- 

monics F,,_, (t?, (p) with ^ = 0, 1 , 2, . . . and m = — £, — {£ — I) ,. . . ,{£ — I) ,£ which come as solutions 
of the Helmoltz equation Af + /(/ + 1)/ = on the unit sphere. 

A convenient basis to describe any function on the sphere is given by these spherical harmonics 
7^^ (7^,9 ) since they form a complete set of orthonormal functions and thus, any square-integrable 
function can be expanded as a linear combination of these. 

fi^,9) = t I c,Y,Ai»,(p), (3) 
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Table 1: Some spaces as limits of the fuzzy sphere. 



where the orthonormalization condition 
allows us to compute the c,^ coefficients as 

We are now ready to define the fuzzy sphere §| of radius R |[T7l[T8l[T9l . It is a non-commutative 
space defined in terms of the N x N matrix operators {xi ,X2,xt,) subject to the relations 

xj+xl+xj = rH, and [xi,xj] = iSijk . --h = ieijk®^, (6) 

VN^ — 1 R 



with = 2R^/ — I and Eijk the totally antisymmetric unit tensor. The operators Xj can be related 
to the angular momentum operators in their irreducible representation of SU (2) of size {2i+ 1) 
with the formula 

= L; = — L;, (7) 

where the relation between the matrix size A'^ and the representation of the angular momentum £ is 
given by A^ = 2£ + 1 . Replacing the equation ([T]) in Q we recover the angular momentum algebra. 

In the table [T] we show some limits of the fuzzy sphere in terms of the matrix size N and the radius 
of the sphere 7?. In that way, the fuzzy sphere contains some other spaces as limits of the matrix size 
and its radius. 

From the algebra of matrices of size N, denoted Mat^, generated by the position operators Xi in 
(O, one can define a Hilbert space, by introducing an inner product. To do that, consider two elements 
of the algebra Mat^ denoted ^ and their scalar product and associated norm are defined by 

Ait Att 

= ^Tr [^^^f] , ||0f = im = -^Tr[(t>'^(t>] , (8) 
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where the normaUzation was chosen so that the unit matrix 1 and the constant function 1 on the sphere 
have the same norm. 

The geometry of the spaces is given through derivation operators. In the case of Sp, the derivations 
correspond to the adjoint action [Li,-] of the angular momentum operators of SU (2). The 
Laplacian is then deduced as 

if2(/)=^-^-0 = [L,-, [!,•,(/)]]. (9) 

Similar to the expansion Q of a function / (t?, 9) on S^, a convenient basis to expand any N xN 
matrix on is the polarization tensor basis. The polarization tensors are denoted by F^,_^ with 
<£ < [N —I) and —l<m< +£, and are defined as the simultaneous eigenvectors of the laplacian 
=Sf^ and axial angular momentum J^y. 

= l{l + 1)F,„„ if3^/,n = mY,„„ (10) 

and we see that .if^ is a cut-off version of — V^. They are normalised to form an orthonormal basis 
of Mat^ 

'-^Tr[YlX,J=5^A,.n (11) 
and transform simply under complex conjugation 

yt=(_i)'«y . (12) 

The expansion of (p in f^,^,, is given by 

N-l +(. 

'/'=I I (13) 

e=Q m=-l 

where the coefficients can be computed by means of the orthonormality condition 

c,^ = _Tr[y>]. (14) 

3 Real scalar field on a fuzzy sphere 

Before introducing the real scalar field theory on the fuzzy sphere, let us look at this theory on an 
ordinary continuum 2— sphere. 

Let be a real scalar field on a sphere with radius R and 0^ potential, the functional action is 
given as 



S[0]= / m 

J§2 



(15) 
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where V,- (/ = 1,2,3) are the usual generators of rotations, r is the mass parameter and A is the 
coupling constant which may depend on the radius of the sphere. 

Second order phase transitions can not appear in finite volume systems, such as the sphere. How- 
ever, it becomes possible in the planar limit, oo. The model on a bidimensional plane, which 
corresponds to the planar limit of the sphere, is defined by the action 



Sm= I d^x 



This model has been widely studied, see for example 12011211 . 

Similarly, the model to study on the fuzzy sphere is aHermitian matrix model which corresponds 
to a real scalar field and is given by the action fSJ |22l 

S[^;N,a,b,c] = Tr \a [L;,0]"^ [L,-,0] +b^^ + c^'^] =Tr[a^\^^^)+b^^ + c^'^] , (16) 

where N is the size of the matrix, b is the real mass parameter, and c is the real, positive, coupling 
constant. Similarly to a commutative sphere, [L,-, •] are the usual rotation generators where L; is the 
angular momentum operator in its irreducible representation of SU (2) with size A'^ = {2i + 1) defined 
by the commutation relations [Li,Lj] = EijkLk- The constant a is a positive number employed to fix 
the units ^ The a term, called kinetic term, contains the information on the geometry of the space by 
means of the Laplacian, while the rest of the action is called potential term and denoted V{^). 
The action ( [T6l) approximates the continuum action ( fT5] ) when 

271 , iTirR^ kXR^ 

These parameters are chosen so that the fuzzy action was normalised so that ^[i] for the unit func- 
tion/matrix be the same on the continuum and fuzzy sphere. 

The absolute minima of this action can be obtained by searching for configurations minimizing 
both the kinetic and potential term separately. The kinetic term is obviously positive and is there- 
fore minimum when _Sf^0 = 0, that is when = ai is proportional to the identity. Replacing this 
constraint in the potential term, we get 

V {ai)=N{ba^ + ca^). (18) 

The necessary conditions to have a minimum are: S' (a) = and S" (a) > 0. If < then we find 
two absolute minima at 



*It is possible to scale ^, b and c to absorb a i.e. fit a to one. The scaling for the field is given by: ^ = Va, leading to a 
scaling for the other parameters of ^ = b/a, and c = c/a^. These changes affect the expectation values by a constant overall 
scaling which has no consequence on such things as phases and phase boundary lines. 
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which have energy S(a) = —Nb^ /Ac, whereas when b > there is only one minimum at a = 0. 
Finally, when b = 0, there is a critical point at a = which is clearly a minimum since S{a) = Nca^ 
and c > 0. 

There are however other local minima to this action which will play an important role in one 
of the phases of this model. They can be located approximately by looking at the minima of the 
potential 111 which are given by U^DU where U is a unitary matrix and D a diagonal matrix with 
diagonal elements ztoo- The absolute minima found above correspond to the particular case when aU 
the diagonal entries of D are identical. 

4 The Metropolis simulation 

We started the simulations by using a standard Metropolis Monte Carlo algorithm |[24l l25l with the 
jackknife estimator for the error |[26l to account for the autocorrelation of the samples.. 

The initial conditions, i.e. the choice of the first configuration in the Markov chain can be of two 
types: Cold initial conditions, which correspond to configurations which are classical minima of the 
action, or Hot initial conditions, which are configurations chosen randomly in the phase space. We 
made sure in our numerical simulations that none of our results depended on the initial conditions, 
whether they were cold or hot. 

In general, when we start the simulation the sequence of samples obtained by Metropolis algo- 
rithm goes through a transient regime where it does not obey the desired statistics yet. This is the 
thermalization process. This is true even in the case of "cold" initial conditions because the classical 
minima may be probabilistically irrelevant when the fluctuations are important. This actually happens 
in one of the phases of our model (the non-uniform ordered phase). 

Tunneling is an essential process in our model as there are multiple classical minima which con- 
tribute significantly to the probability distribution of the field. Typically, tunneling is exponentially 
suppressed by the energy barrier separating the classical minima. It can therefore be difficult to 
account for in the Monte-Carlo algorithm. To improve the probability of tunneling, we have tried 
various sampling methods. 

The two simplest ways of sampling the phase space are to either make a big change on the matrix 
as a whole or to perturb its entries one by one. The first method allows for big changes and helps 
tunneling but usually yields unfavored, high energy, test configurations which aie rejected and in- 
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crease the autocorrelation between configurations. On the other hand, the latter is good at exploring 
the phase space locally, but has a low chance of tunneling. Even alternating the two methods to enjoy 
both their advantages is not sufficient to produce the tunneling necessary in the model studied. 

As we already discussed at the end of section [3l the classical minima of the action are located at 
=b ^J —hjlc. Thus, the interval where we must vary the real and imaginary part of every entry of the 
matrix during the sampling should be about / = \—l\J —b/2c;2\/ —b/2c\. In practice, we have found 
empirically that we need an interval of variation of the field between 2.3 and 2.6 times bigger. 

When we use an interval less than 2.37, the trace effective probability density distributions of 
the matrix will not reproduce the results obtained via direct integration for N = 2. In general, this 
effect also appears for any matrix size N. The upper bound does not affect the results so much as the 
auto-correlation of samples (more configurations are rejected by the metropolis algorithm) and thus 
the speed of convergence of the code. We have found that 2.6/ is the optimum upper bound to balance 
speed and precision. 

A more sophisticated method we have successfully implemented is the annealing method ||25]| . 
The idea is to produce favored decorrelated test configurations by introducing a temperature-like pa- 
rameter j8 in the probability distribution exp(— jSS'). The Metropolis sampling is done normally with 
j8 = 1. During that sampling, the field is typically trapped around one of the classical minima. Peri- 
odically, the Metropolis sampling is interrupted, and this parameter is lowered (i.e. the temperature 
is increased) which smoothes out the action and allows the field to move more freely between the 
classical minima. Then j8 is raised back to one (i.e. the temperature is lowered back) trapping the 
field around a classical minimum which is hopefully decorrelated from the previous one. 

The annealing method thus increases the probability of tunneling between minima of the action 
and decreases the autocorrelation between configurations. It also increases the computation time, but 
the gain in efficiency is largely dominant, making this method very useful and reliable for simulations 
with larger matrices. 

The computation time can still be too large, making the simulation impossible to run in practice. 
We have developed a method where the real time of computation decreases dramatically which we 
present in the appendix. 

5 Limiting models 

In this section we present the lowest dimensional model which can be integrated directly and the pure 
potential model which can be solved analytically. They will both be useful to test the validity of our 
Monte-Carlo simulation in some particular limits. 
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5.1 Lowest dimensional model 



It is quite useful to investigate the lowest dimensional model {N = 2) as it has only two independent 
parameters and can therefore be well understood, and integrated directly. This provides a good in- 
dependent computation to test our Monte-Carlo code against. Furthermore, it happens that even this 
low dimensional case shows all the features of the large limit! 

In the simplest case N = 2, the action (fT6l ) can be simplified by expanding the field in terms of an 
orthonormal basis {l,a^}, where 1 is the 2 x 2 identity matrix and a<; are the three Pauli matrices. 
The expansion is 

= al + • "a (20) 

where the coefficients a and p^. are in M. Then, writing down the action ( fT6l ) in terms of this new set 
of variables, we get 

S = 4ap^ + 2b (a^ + p^) + 2c (a^ + 6a^p^ + p^) (21) 

where p is the norm of p^. 

The action (|2T] ) depends only on the modulus of the vector 'p. This property allows us to inte- 
grate out the degrees of freedom associated with the rotational symmetry of the vector 'p, which is 
the expression of the general SU (2) symmetry of the action in two dimensions. The corresponding 
effective probability density distribution is given by 

P^s[a,p] = ^p^e-sM = ie-Seff[«,p] ^ z = J dadpp^e"^ = J dadpe^^^", (22) 

where S^ff = 5 — In p^ is the associated effective action. 

This simple example depends on two variables only, which makes it possible to integrate numer- 
ically without a lot of effort for any set of parameters {a,b,c} via the trapezoidal rule or any other 
algorithm ||27l . to get the expectation values. In this sense, we can solve directly the model for any 
set of parameters making it possible to test our Monte Carlo code. Any graph in this section will not 
include error bars because, with direct integration, they are negligible. 

Because of the SU (2) symmetry of the theory, the expectation values of (|a|) and (p) give us the 
whole information about the average configuration (0). We can see their behavior in the figure [T] 
computed from a direct integration with N = 2, a = \, and for a typical value of c = 50, as a function 
of the remaining parameter b of the model (scaled to bc^^^^). We can see three distinct phases: 

1. Disordered phase: the expectation values of |a| and p are close to zero, roughly in the interval 

(0,+oo). 

2. Uniform ordered phase: the most important contribution to the configuration is given by the 
expectation value of a, roughly in the interval {—°°, —24) 



Q 
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Figure 1 : Expectation value of | a | and p obtained from direct numerical integration. 

3. Non-uniform ordered phase: the system is ordered but the main contribution to the configuration 
is given by the expectation value of p, in between the previous two phases. 

The disordered phase has averages of |a| and p, in the expansion ( [201 ). at approximately and 
the typical configuration is thus distributed around zero. This phase is analogous to the paramagnetic 
phase in ferromagnetic materials. Following the analogy, if we take the parameter as a ''temperature 
parameter", the thermal fluctuations do not allow any kind of ordering in the material when the 
temperature is bigger than some critical value. The thermal fluctuations are getting stronger and 
stronger when the temperature is increased. 

The uniform ordered phase is characterized by the fact that the most important contribution to 
the configurations is given by the coefficient a, in the expansion (l20l ). The expectation value of 
p is negligible with respect to the expectation value of |a|. This means that the configuration is 
approximately proportional to the identity matrix. This phase is analog to the magnetic phase, in 
ferromagnetic materials. 

In the third and last phase (the non-uniform ordered phase), both |a| and p contribute to the 
configuration but in this region of the parameter space, p is more important than |a|. This phase 
has ordering i.e. the field has non-zero expectation value but this ordering is not an analog of any 
ferromagnetic ordering. It was argued in ifTTl \V2\ that in this phase, the eigenvalues of the matrix 
has two cuts located at the two minima of the action itao given in ([T9l ). whereas IH speculated that 
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N=2 a=1 , b=-257.5, c=50.0 30.03.05 



Trace/(-bN''/2c) " 



(a) Probability distribution of a. 



N=2 a=1, b=-257.5, c=50.0 30.03.05 



p/(-b/2c)^'^ 

(b) Probability distribution of p . 



Figure 2: Comparison of the unnormalised probability density distributions of some observables for 
A'^ = 2. In most cases, the two curves cannot be distinguished. 



the eigenvalues of would be split equally between positive and negative eigenvalues. For N = 2, it 
means trivially that (/> = OoCJa up to a free SU (2) rotation, and thus one would expect < |a| >^ ao 
and < p >~ Oo = y/-b/c^-^ /{4c)^-^^ . This is indeed true in figure [B as for -20 < b/^/c < -5, 
< p > does curve like y/—b/c^-^/ 200*''^^ ~ ^/—b/c^-^/4 and is much bigger than < |a| >. 

As stated earlier, we can also use the results from this alternate method to validate the Monte- 
Carlo code for N = 2. The figure |2] shows the unnormalised effective probability density distributions 
of the quantities a and p. In that case, we can compare directly both effective probability density 
distributions. The excellent agreement shows that the statistical error bars are negligible, but the most 
important thing is that the program has sampled the phase space properly. We have already checked 
many points in the parameter space and we have always obtained identical results up to error bars. 
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At the moment, we have shown the convergence of our simulation in the lowest dimensional 
model, but our goal is of course focused on simulating the model using bigger matrix sizes to extrapo- 
late to the continuum limit {N +°°). Still we will see that the N = 2 results are already remarkably 
good approximations of the large limit. 

5.2 Pure potential model 

The pure potential model interests us for two reasons: it can be solved analytically and gives a good 
approximation for the transition curve between the disordered and non-uniform ordered phases dis- 
cussed previously HI. It comes from setting a = in the action ([T6l ). only keeping what we called 
the potential term. This approximation is increasingly accurate as the transition is tracked to larger 
couplings far from the triple point. 

This model {a = 0,N ^ o°) has been solved by many authors ||28ll29ll . In term of their solution 
we can get an expression of the specific heat and other thermodynamics quantities which are a good 
reference to compare to the numerical results and the convergence of the algorithm when we increase 

The specific heat in this approximation has the form 

^ '^"^ (23) 

i+^_^(2r2-3)v^ f>-l 

where r = b/\bc\ with be = —2^/Nc is the critical mass. From equation (1231 ) the phase transition is 
a third order transition because the first derivative of the specific heat has a finite discontinuity in 
b = -l. 

Another way to detect the phase transition is to look at the probability distribution of the field 
eigenvalues. In the disordered phase, they are confined into a single connected region centered around 
zero, whereas in the non-uniform ordered phase, they are split into two disconnected regions centered 
respectively around ±y^—b/2c corresponding to the minima of the polynomial potential. Due to this 
characteristic behaviour, we also refer to this as a "one cut-two cut" transition. We will also use this 
terminology for the disorder/non-uniform transition since the work of Panero ifTTl [T2l shows that the 
transition in the fuzzy sphere model, where a y^O, also have this characteristic behaviour. 

The phase boundary for this model is given by 

b = -2VNc (24) 
and is included in the phase diagram shown in figure [6] at the end of this article. 
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Figure 3: Specific heat for a = 1 and c = 40. 

6 Observables and Simulations 

For the model under study, the number of degrees of freedom is A^^, which corresponds to the number 
of independent real entries in a Hermitian matrix. Thus, the thermodynamic limit we are interested 
in corresponds to matrices of infinite size. The standard procedure, to take the thermodynamic limit, 
is to define a scaling of the parameters of the model such that the relevant observables collapse in a 
phase diagram independent of the matrix size. If the phase diagram collapses in a reasonable way, 
then we can straightforwardly extrapolate it to the thermodynamic limit. 

The specific heat is a measure of the dispersion of the energy. It is sensitive to the phase transitions 
which register as peaks in it. We therefore use it as the order parameter. Typically, it will present one 
or, more often, two peaks as we show in figure [3] for {a = 1, c = 40} and various matrix sizes. The 
very obvious peak is located around bN^^/'^ = —6 ±0.4, the other one, almost imperceptible is around 
_ _2 Pqj- jj^g biggest matrix size investigated N = 10, simulations for a curve as the one in 
figure [3] took about a day. The error bars provided by the jackknife algorithm were omitted as they 
are quite small and would only crowd the figure more. 

Other observables, such as < Tr [0^] > and < |Tr [0] | > which were used as order parameters in 
ll6l . their susceptibilities, and the internal energy < 5 > have also been collected but are not shown 
here. They will be used in section |7] to identify the phases though. 



It is an important remark that the transition, from the non-uniform ordered phase to the uniform 
ordered phase, presents a very high and wide peak in the susceptibilities which can subsume and hide 
the smaller one when near the triple point, making it impossible to determine its exact position. As a 
result, the data points of this transition curve in the phase diagram could not be found near the triple 
point. However, Panero lITTl . by looking at the eigenvalue distribution of for c/aN'^ = 1/2 provides 
an additional point on this curve very near the triple point. 

In the figure [3l the scaling for b, given by bN^^/'^, which ahgns the peaks (and thus the location 
of the phase boundary) for different matrix sizes has already been included. It is remarkable that with 
this scaling, the N = 2 curve has the same quahtative behavior as the A'^ = 10 curve, as announced 
previously, but the peak in figure [3l is already a reasonable approximation of the large limit peak 
found for N = 10. 

This analysis was repeated for a wide range of the parameter c and for matrix sizes < 10. The 
collected results, the interpretation of the phases and the collapsed phase diagram will be presented 
in section |7J 

7 Results 

In this section, we wiU present the collapsed phase diagram as well as an analysis of the three phases 
observed. 

In the plots |4] and [51 we can see different profiles of the probability density distributions as a 
function of the mass parameter b for Tr [0] and p = ^y\c^ + |c,op + |cj| p which gives the power 
in the 1 = 1 angular momentum mode in (fT3] ). with {a = I, c = 40, A'^ = 4}. There, we can appreciate 
the three different phases of the model. 

In the uniform order phase, for b negative enough, the trace is distributed around two symmetric 
values centered on ibao respectively, ^, and p is distributed close to zero, i.e. it gives no contribution 
to the typical configuration. In this phase, is approximately proportional to the unit matrix and the 
rotational symmetry is thus preserved. 

The non-uniform ordered phase, for intermediate values of b, has the peculiarity that the most 
exterior peaks of the probability of the trace, which correspond to the absolute minimum of the action 
±0£o and thus to the field in the uniform order phase, are smaller than the new peaks which arise 
between them. Furthermore, the most probable value of p is not close to zero. In this phase, the 
power of the configuration is thus in higher angular momentum modes (as defined in the expansion 

was defined in ( fT9l l as the location of the absolute minimum of the action. 
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N=4 , a= 1, c= 40 




Figure 4: Probability distribution of the trace of as a function of b 



N=4 , a= 1, C=40 



Probability Distribution 




Figure 5: Probability distribution of p as a function of b. 
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Figure 6: Full collapsed phase diagram. 



(fT3l ) in polarization tensors) and the rotational symmetry has been spontaneously broken. 

The last curve is representative of the disordered phase. In this phase the configurations (both the 
trace and p) are spread over a long interval but very close to zero restoring the rotational symmetry. 

A phase diagram is a map that contains the thermodynamics or physical properties of a given 
system. This imphes that, to construct a phase diagram, we need quantities in the thermodynamic 
limit. As explained in section |6l this is done by finding a scaling in the bare parameters of the model, 
b and c here, to make it independent of the number of degrees of freedom A^. 

We had already found in section |6] that the scaling necessary to make the diagram independent of 
N was 5 for the mass parameter b. Repeating the simulations for various values of c and plotting 
the phase boundaries found for all values of N simulated, we found a scaling in N^-^ for the coupling 
constant c. This scaling is the same for all the coexistence curves which guarantees a consistent 
N ^ o° limit. We can then define scale-free parameters 

— b _ c 

b = — tm ) c = . (25) 

Remember that for all the simulations and results in this paper, we have set a = 1 . 

These results are presented in the figure [6] which shows the phase diagram for the 0^ model on 
a fuzzy sphere. The three phases we identified above are delimited by the coexistence curves which 
meet at a triple point. These coexistence curves can be fitted to get an algebraic expression for each 
one of them using the scale-free parameters introduced above in ( [25] ) 
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As mentioned in section [6] we could not access the Disorder/non-uniform order phase boundary 
near the triple point. However since the curve is consistent with a straight line, we can extrapolate it 
to the triple point without any difficulty. We find: 

Disorder/non-uniform order: c = 2.29{—b) —4.74. (26) 

As expected, for large c, this curve is well approximated by the one obtained for the pure potential 
model derived in |5.2[ given in scale-free parameters by 

c={-W/^, (27) 

and drawn with a dashed hne in the phase diagram, figure [6l 

We did not focus on the disorder/uniform order boundary line in this paper since it has already 
been studied in detail in 161. It was found there to be a straight line going through the origin. Con- 
verting its equation to our scale-free parameters through (|25] ). we get 

Disorder/uniform order: c = 0.23(— Zj). (28) 



Finally, the uniform-non-uniform order phase boundary line which was studied in detail in this 
paper is approximately straight with equation 

Uniform - non-uniform order: c = 0.2{—b) + 0.07 (29) 

which just prolongs the disorder/uniform order line, up to error bars. 

These three coexistence curves, ( I26I28I291 ). intersect at a triple point given by 

{bT,CT) = (-2.3,0.52). (30) 

These values are consistent with the data presented in ifTTI . In fact, figures (11-30) there correspond 
precisely to c = 0.5 ~ cj, and by identifying the point where the eigenvalue density undergoes the 
one cut-two cut transition described in section ls!2l one finds that his data gives br ^ 2.3 consistently 
for A^= 15,17,19,21,23. 

If instead one takes the asymptotic form of the disorder non-uniform order transition line given by 
the one cut-twocut transition ( |27] ) instead of (|26l ), and finds its intersection with the disorder-uniform 
order transition curve (|28] ). the triple point occurs at 

{b^T,c'r) = (-0.92,0.21). (31) 

We conclude from this that the effect of the kinetic term is to move the triple point to larger values 
along the line governing the disorder/uniform order transition. 
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8 Conclusions 



In the main part of the paper, we presented the results for the numerical simulation by means of an 
optimized Metropolis algorithm for the 0^ matrix model. In the appendix we develop the metropolis 
algorithm which makes more efficient the simulation of matrix models. In particular, we argue that the 
algorithm proposed presents considerable advantages with respect to the usual Metropolis algorithm 
in the simulation of matrix models |f30|. The reduction in the processing time for both non-uniform 
ordered and uniform ordered phases will be more evident for large matrices and, of course, when we 
are far away from the coexistence curves due to the fact that the minima of the potential are more 
separated. A different approach was used with equal success in |fTTT|. 

Figure |6] shows the phase diagram for the model given by ([T6l ) and refines the phase diagram 
which was incomplete in Q. The data have been collapsed using the scaling form shown on the axis 
and defined in (1251 ). It is consistent with the scaling of the exact solution of the pure matrix model 
which only fixed the quotient of the two scalings. One of the important features of the diagram is that 
all three coexistence lines can be collapsed simultaneously. This did not have to happen and in fact 
the corresponding lines do not all collapse together for a related three dimensional model IIBTTI . where 
the spacetime is taken to be a fuzzy sphere direct producted with a temporal direction. 

This diagram contains the information about three different phases, the well known disordered 
and uniform ordered phase, and a new phase, the non-uniformly ordered phase (where the S0{3) 
spatial symmetry of the round sphere is spontaneously broken), as well as the scaling of the model, 
and the coexistence curves. We could even estimate the coordinates of the triple point which is the 
point where the three phases coexist in equilibrium. The coordinates of this triple point are consistent 
with the independent simulation ifTTTl . 

Another article llT3l finds different results, including an extra phase and no scaling. An obvious 
reason may be that in the phase diagram they show for N = 25, our scale-free parameter —b has a very 
small range in our scale-free parameters of [0;0.13], meaning it only shows a tiny corner of our phase 
diagram of figure [6l Furthermore, they use the probability distribution of (pn (denoted On there) as 
an observable to detect the transition between the two ordered phases. First, this does not seem to be a 
physically meaningful observable, especially given the SU (2) symmetry of the model. Furthermore, 
they obtain a curve somewhat similar to the one cut curve for the eigenvalues of 0, but they locate the 
transition when this profile gets deformed with a dip at zero, instead of when it switches to two cut (if 
it ever does). This boundary line is absent in our phase diagram and in previous ones ||6lllI]|T2l, and 
we find no evidence for such a transition or a new phase in this region of the phase diagram. As for the 
lack of scaling for the other phase boundaries, it is difficult to decide the cause, but it is disquieting 
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that their simulations sometimes depend on the initial condition, such as when they find hysteresis. 

In the large limit the model with a = has a third order phase transition between disordered 
and non-uniform ordered phases ||28ll29l . The disordered phase is described by a single connected 
eigenvalue distribution called a "one cut phase" distribution, whereas the ordered phase is described 
by an eigenvalue distribution split into two disconnected distributions centered on opposite values 
and called a "two cut phase" distribution. The transition occurs when the two cuts merge to become 
a single cut for c = /AN. Figure |7] confirms numerically the convergence of the disordered/non- 




Figure 7: Plot of the specific heat at the disordered/non-uniform ordered transition for increasing and 
its ^ oo limit, the exact pure potential model, given by Eq. (|23l ). 

uniform ordered transition towards this exact critical line of the pure potential model as the coupling 
is increased. The simulations of Panero ifTTl [T2]| confirm that this transition for the full model is 
indeed a one cut-two cut transition though the eigenvalue distribution now has a richer structure. 

We expect that the existence of the cut transition of matrix models and of a triple point is a generic 
feature of fuzzy scalar field models, since all such models should reduce to a pure matrix model when 
the kinetic term becomes subdominant. This means that fuzzy scalar field models should generically 
have an exotic phase with spontaneously broken spacetime symmetry. 

Numerically, it is not difficult to find the coexistence curve between the uniform ordered and 
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disordered phases which exist for low values of c. On the other hand, the coexistence curve between 
the two ordered phases is difficult to evaluate because it involves a jump in the field configuration and 
tunneling over a wide potential barrier. 

In the current model the triple point is estimated to be located at 



This is obtained by extrapolating the three coexistence lines till they meet. Surprisingly good agree- 
ment with this result was obtained in fZl by performing perturbation theory in the kinetic term, i.e. 
by expanding in the parameter a to second order. It is not however, totally clear that the triple point 
identified there coincides with the one here as a different scaling of the parameters was necessary, but 
the salient features are the same. 

The position (|3TI ) where the curve (l27l) intersects the disordered /uniform order transition line 
suggests that the effect of the kinetic term is to push the transition, and hence the triple point, to larger 
negative values of b. This is a positive feature since it indicates that adding a higher derivative term 
to the model will allow one to tune the triple point to large coupling. The conjecture is that this will 
be sufficient to eliminate the UV/IR mixing problems ||22]| and recover the commutative theory with 
the correct fluctuations |[23l . 

It still remains to be seen what thermodynamic limits can be drawn from the phase diagram and 
the scaUngs in each of the limits of the fuzzy sphere introduced in table [T] the ordinary sphere, the 
ordinary plane, and the Moyal plane. 

To that end, we want to reexpress the positions of the coexistence curves (I26I28I29I ) which depend 
on a = \,b,c and A'^ as a function of parameters well defined in the thermodynamic limit. These are 
the radius of the sphere R and the non-commutative parameter = 2R^ / VN'^ — 1 appearing in Q 
and in the list of possible hmits of the fuzzy sphere of Table [T] and r and A appearing in the action 



(^r,cr) = (-2.3,0.52). 



(32) 



(fT6TT7] ). 

Using the scalings of (1251 ) and N ~ 2R^/&, we find 

- r03/2 _ A0 



(33) 



2V2R' 4871' 

which can now be replaced in the algebraic fits for the coexistence curves to get 




(34) 



Disorder/uniform order: 




(35) 




(36) 



(37) 



Since the phase boundary lines and the triple point all scale in the same way, it is not surprising to 
find that, out of the four physical quantities available, only two are independent: r&^^^/R and A0. 
As a result there are not enough physical parameters to fix the limiting procedure. For instance, if the 
limiting space, represented by R and 0, is fixed, one can still scale the field model parameters r and 
A freely. 



A Optimized algorithm for matrix models 

We now present an improved Monte-Carlo scheme we used to speed up our simulations. 

The probability transition function (denoted by PTF for short) of a Monte-Carlo algorithm W/ 
from an initial state / with probability P, to a final state / with probability Pf, must satisfy the detailed 
balance equation 

PiWf,i = WijPf. (38) 

The Metropolis PTF 

r Pfi 

(39) 



Pf 
1,— 
' P 



Wfi = min 

is the best known example of one such, but another introduced in | 32] is given by 

(40) 



Wf = w, min 

fj J.I 



1 

' P' ^/.. 



which is a generalization of the Metropolis PTF when w^. 7^ w.^.. In our case, we have selected the 
further generalization 

Pf' 



Wj, . = mm 



Pi 



(41) 



which is equivalent to the Metropolis PTF using a different probability distribution p yet to be defined. 

The Boltzmann probability density distribution P (a) to find a configuration in the volume {x,x + dx) 
is defined by 

P(x) = i^'-^H, (42) 

where S [x] is the action or energy, and Z is the partition function which contains the whole relevant 
information of the system. In general, it is not possible to obtain an exact expression for the partition 
function analytically or numerically. The Monte-Carlo algorithm via the Metropolis PTF (|38] ) is so 
important because it does not depend on Z. 

With AS = S[xf]— S [xj], putting (|42] ) in (|39l ) gives us the PTF in terms of the Boltzmann proba- 
biUty density distribution (denoted PDDfor short), that is 

Wf J = mm[l,e-^^]. (43) 
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AS<0 


AS>0 


<0 
>0 


Wf.=mm[l,e-^^+^'] 
wf. = e-^' 


f-.i 

l^B ^min [e-^',e-^^] 



Table 2: Possible cases in the evaluation of the Probability Transition Function Wf.. 

In the same way as (l42l) . we can associate a kind of energy s (y) to the probability distribution p 
introduced in (|4TI) 

p{y) = -e-'^l (44) 
z 

Now, the new Metropolis-Boghosian PTF which comes from the equations (l40l ). (|4TI ). (l42l ) and (l44l ) 
is given by 

Wl = min [1,^'-^-^] min = min [\,e-^\e-^' ^e-^'^"-'] , (45) 

where Ai' = s [xf] — s [x,] is the equivalent of the AS defined above. All the possible cases for the PTF 
( |45l ) are presented in the table |2l 

We can ask why we might need the PTF (l45l) when we have a simpler function (1431 ) already? 
When the evaluation of AS is quite simple, for instance for the Ising model, this methodology is 
counterproductive because more exponential functions must be evaluated. On the contrary, when the 
evaluation of AS is computationally very expensive, as the matrix models are, the PTF (l45l) avoids 
the evaluation of very improbable changes in the configurations due to the implementation of the filter 
As, which reduces the processing time. 

Numerically, we do not want to evaluate both AS and A^. If 5 is a good approximation of the 
effective potential created by S but simpler to evaluate then, we can use the Metropolis algorithm 
with the action s to refuse or accept the new configurations before the evaluation of AS (which is 
complicated to evaluate and only will take machine time). 

In this section, we propose a variant calculation of PTF (1451 ) shown in table |3l where we avoid the 
evaluation of AS when the previous evaluation of Metropolis algorithm with As refuses the attempt to 
change the configuration. 

A.l Relative error 

The new PTF, denoted by Wj., can not satisfy the detailed balance equation ( [39] ). This fact introduces 
deviations in the probabilities, in exchange for a computational time reduction, since A^ will be chosen 
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i 


A5<0 


A5>0 


< 
> 


Wj.= mm [1,6-^+'^] 


Wj. = e-^ 



Table 3: Proposition for a faster Probability Transition Function. 





A5<0 


A5>0 


A5<0 












J ifAS<A5 


A5>0 











1 l_e-A5+A.T ifA5>A5 



Table 4: Relative error (1461) between the probability transition functions, W^. and .. 



simple to calculate. Anyway, we will make sure to keep under control the error introduced by this 
breaking of the detailed balance equation. 



In accordance to this, the relative error between Wf and W, is defined as 



err 



wf 



(46) 



and its values are shown in table ID As we can see in that table, only the case A5 > As > presents a 
relative error different from zero. This error goes to zero when A5 > As and it goes to one when A5 
As. Similarly, when A^ » 1 we can almost take for granted the rejection of the new configuration by 
Metropolis. Thus, the introduction of the PTF Wj. is very convenient to estimate the PTF given by 
(|43] | breaking the detailed balanced equation, where we only expect a tolerably small deformation in 
the averages (specifically in regions with low probability) with respect to the averages obtained from 
the PTF's (|43]l and (|45]l. 

In our experience, we can obtain a better approximation when we replace As {AS)' only in 
the case A5 > As > 0. The prime indicate the difference of energy obtained by Metropolis one step 
before under the condition A5 > As > 0. In that way, the algorithm has a kind of "auto-regulation" 
which reduces the deviation of the averages with respect to the PTF which obeys the detailed balance 
equation. With this auto-regulation we only update the energy reference level. 
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Figure 8: The FDD's obtained by different methods for N = 2. 

A.2 Application to the fuzzy scalar field model 

We can now adapt this scheme to the matrix model (fT6l ) we are considering in this article. A finite 
variation of the action, from a configuration to a configuration (p + (p is defined by AS = S [<p + — 
S[(j)], where (j) must be a Hermitian matrix. The evaluation of (A^)^^, for the (j)^ matrix model for a 
single entry (/x,v), involves the evaluation of a cubic polynomial in the matrix. This is a highly 
non-local function which causes the main slow-down in the code. 

Figure [8] shows a comparison for the same set of internal parameters (Monte Carlo time, ther- 
malization time, decorrelation time, etc.), at a collapsed point {b,c) = (—2^/^,1) from the phase 
diagram -Figure IH- corresponding to the non-uniform ordered phase, between the results obtained by 
direct integration (explained in Section ISTI ). and the Monte Carlo simulations via either of the three 
probability transition functions presented above: Metropolis, Metropoiis-Boghosian, and Metropolis- 
Boghosian-Fergar (the one in Table O. The small deviations (noise) of the Monte-Carlo simulations 
with respect to the direct integration are a normal effect of the Monte Carlo simulations and can be 
reduced by increasing the number of samples produced by the code. 

Going back to the choice of the filter s introduced in (l44l ). we want a function that incorporates 
the main features of S but is easier to evaluate. 

We can notice that the action (fT6] ) for a fixed entry (/i, v) of the matrix ^ correspond to a quartic 
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polynomial in the entry ^^y, thus 

where the coefficients A, B and C depend on the rest of the entries in the matrix. In that sense, we 
propose 

V=C>Mv+fi>M'v+^^ (47) 

where A', B' and C' are constant coefficients. Thus, s^^ goes to S^^ when {A',B',C'} {A,B,C}. We 
will obtain a better concordance between both AS' and by choosing a set of parameters {A',B',C'} 
close to the non-primed parameters. As a first approximation, we took 



{A',B\C'} 



{0,b,c} if /I = V 



{0,0, c] if Ai/v 



where b and c are respectively the bare mass and interaction parameters of the model ([T6] ). For 
simplicity we have fixed A' = but in general, we can use any other real number and it will not affect 
As. This set of primed parameters for s, was chosen to contain the most basic information of the full 
model S. 

At the end of section [3l we have shown that the action (fT6l) has two symmetric minima with 

h_ 
'2c- 

a with =±Y^2c' 

(pi = \,2, . . . ,N). K simple function of with the same set of minima has been given in the equation 
(l47l l with the parameters (A' = 0,B' = b,C' = c). 

For non-diagonal entries, we have observed that their probability distribution is around zero. Thus, 
it is enough to consider the function ( |47] ) with the parameters (A' = 0,B' = 0,C' = c). 



respect to the trace. Those minima are located in Tr [(/>] = -^N \J — ^. Thus, we can consider that 
every single diagonal entry in the matrix contributes to the trace minima with = ±1/— j-, where 



A.3 The algorithm 

The algorithm is basically the same as the usual Metropolis algorithm although with some adaptations 
to the current setting. In the figure |9] we show the flow chart for the implementation of the new method 
that we have proposed. 

Internal variables in the flow chart |9] 

• X and x': random numbers uniformly dis- • Metropolis(A/): indicates the Metropolis 
tributed in the open interval (0, 1). algorithm using the difference of energy 



OS 




Figure 9: Flow chart of the Optimized Monte-Carlo method. 



A/. • p: ratio of implementation of the new 

. j8: represents the difference of the energy "^^^^^^ respect to the standard one. 

from a Monte-Carlo step which had been • (T): part of the subroutine where we avoid 

evaluated before. to evaluate AS'. 

In the flow chart presented above, we have emphasized with (T), the step in the simulation where 
we avoid the evaluation of A5. This would seem to be insufficient to reduce in a significant way 
the processing time but it is not true at all. The number of times that the algorithm passes trough 
(T^ divided by the total number of times that the Modified Metropolis Algorithm (MMA for short) 
has been used, will be an estimation of the efficiency of the new method with respect to the Usual 
Metropolis Algorithm (UMA for short). 

Let us first define an efficiency parameter for our MMA to compare it to the UMA. If is the 
Monte Carlo time to run the simulation for N xN Hermitian matrices under the model ([T6l) . and T the 
number of times that the algorithm passes trough the new feature (T^, then 

eff=^, (48) 

defines the efficiency of the modified algorithm. In particular, if ?f^„ is the time for a run with a UMA, 
which is without our routine (T) ^nd t the time with it then, we have f ~ ( 1 — eff ) fj^,, . 

Remember that the figure |9] only represents one attempt to switch one entry and, if (5 has been 
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updated then it must be saved for the next one. 

The modification of the Metropolis algorithm presented in this chapter allows us to simulate ma- 
trix models with a decrease of the calculation time with respect to the usual method. The explicit 
breaking of the detailed balance equation by our proposition involves a systematic error which we 
can keep under control at any time. 



A.4 The optimized Metropolis method 

As a second part, we present the results obtained with the optimized Metropohs method. As explained 
in section IAT2I this method was successfully tested in the lowest dimensional N = 2 case. 

As we saw in the figure |9l the Modified Metropolis Algorithm (or MMA) had to evaluate three 
exponential functions compared with the Usual Metropolis Algorithm (or UMA) where we only have 
to evaluate one. Thus, when the efficiency eff defined in ( |48] ) is too small, it could be better to use the 
UMA. This happens for instance, in the disordered phase: the difference in processing time between 
UMA and MMA is not appreciable^. Even worse, the processing time in MMA could be a little 
bigger in that phase. 

It is not the same for the other two phases where tunneling plays an important role. There the 
efficiency eff goes to one and the MMA is greatly more efficient"*. 

To keep under control the relative error when eff is close to one, we have to adjust the p ratio. 
Thus, p « means a fast run, but could present a considerable relative error. At the other extreme, 
when ~ 1, the run will be slow but the relative error will be very small. We have to look for a 
balance between accuracy and speed. We have set p between 0.55 and 0.70 but it is also possible to 
set it dynamically. 

As an example, in the figure \10\ we show the behavior of the processing time per configuration 
with respect to the matrix size obtained by means of the usual Metropolis algorithm for some given 
processor^ when a= \, b = —4 and c = 0.10 which corresponds to the uniform ordered phase where 
we expect some gain. And indeed, the best fit for this curve time = (1.49 ±0.02 x \0^^)N^ + (3.18 ± 
0.02 X 10"^)A^'*s grows like A^'^. 

^They have approximately the same velocity of processing because a large percent of attempts will be in the range of 
fluctuations of s. 

^In these phases, the new method is faster than the old one because a large percentage of attempts could be out of the range 
of fluctuations of s, thus avoiding the evaluation of A5. 

^In this example we have used a Mobile Intel(R) Pentium(R) III CPU - M 800MHz, 369.10Mb RAM. On gcc-4.0.2 2005- 
10-01, Ubuntu, kernel 2.6.12-10-386. Kubuntu 5.10 Breezy Badger. 
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Figure 10: Processing time per configuration for the UMA. 



Starting from a random configuration, it can be thermalized or decorrelated using the MBF 
method described by Table [3j then we can use the usual Metropolis algorithm to evaluate the prob- 
ability of transition between the old configuration and this new sample obtained from MBF. Doing 
this, we save processing time and, at the same time, we do not introduce any systematic error because 
the usual Metropolis algorithm will reject or accept the new configuration which only contains the 
statistical error. 



References 

[1] J. Hoppe, MIT Ph.D. thesis (1982); J. Hoppe, Elem. Part. Res. J. (Koyoto) 80 145 (1989). 
H. Grosse and P. Presnajder, Lett.Math.Phys. 33, 171 (1995); H. Grosse, C. Khmcfk and 
P. Presnajder, Commun.Math.Phys. 178,507 (1996); 185, 155 (1997); H. Grosse and P Pres- 
najder, Lett.Math.Phys. 46, 61 (1998) and ESI preprint, (1999); H. Grosse, C. Khmcfk, and 
P. Presnajder, hep-th/9602115 and Commun.Math.Phys. 180, 429 (1996); H. Grosse, 
C. Klimcik, and P. Presnajder, in Les Houches Summer School on Theoretical Physics, 1995, 
hep-th/9603071. S. Baez, A. P. Balachandran, S. Vaidya and B. Ydri Commun.Math.Phys. 
208, 787-798 (2000); . P Balachandran, T.R.Govindarajan and B. Ydri, Mod.Phys.Lett. A15 1279 



9,8 



(2000); G.Alexanian, A.P.Balachandran, G.Immirzi and B.Ydri J.Geom.Phys. 42 28-53 (2002); 
Badis Ydrihep-th/0110 00 6. 

[2] J. Medina, W. Bietenholz, F. Hofheinz and D. O'Connor, PoS LAT2005, 263 (2006) [arXiv:hep- 
lat/0509162]. 

[3] H. Steinacker, JHEP 04 (2004) 077 [arXiv:hep-th/0501174]; H. Steinacker, [arXiv:hep- 
th/0511076vl]. 

[4] J. Ambjom and S. Caterall, Phys. Lett. B 549 (2002) 253 [arXiv:hep-la1/0209106]. 

[5] W. Bietenholz, F. Hofheinz, and J. Nishimura, JHEP 0406 (2004) 042 [arXiv:hep-th/0404020]. 

[6] X. Martin, JHEP 0404, 077 (2004) [arXiv:hep-th/0402230]. 

[7] D. O'Connor and C. Saemann, JHEP 0708 (2007) 066 [arXiv:0706.2493 [hep-th]]. 

[8] D. O'Connor and C. Saemann, arXiv:0709.0387 [hep-th]. 

[9] H. Steinacker, JHEP 0503 (2005) 075 [arXiv:hep-th/0501174]. 

[10] F. Garcia Flores, D. O'Connor and X. Martin, PoS LAT2005 (2006) 262 [arXiv:hep- 
lat/0601012]. 

[11] Marco Panero. JHEP 0705 (2007) 082 [arXiv:0608202[hep-th]]. 

[12] M. Panero, SIGMA 2 (2006), 081 [arXiv:hep-th/0609205]. 

[13] C. R. Das, S. Digal and T. R. Govindarajan, arXiv:0706.0695 [hep-th]. 

[14] A. P. Balachandran, B. P. Dolan, J. H. Lee, X. Martin and D. O'Connor, J. Geom. Phys. 43, 184 
(2002) [arXiv:hep-th/0107099]. 

[15] B.P Dolan and D. O'Connor, JHEP 0310 (2003) 060 [arXiv:hep-th/0306231]. 

[16] Jurgen Jost, "Riemannian Geometry and Geometric Analysis". Spring-Verlag 3rd edition. 
Berhn, 2002; Charles W. Misner, Kip S. Thome, John Archibald Wheeler, "Gravitation". W.H. 
Freeman. New York, 1970. 

[17] G. Landi, arXiv:hep-th/9701078. 

[18] M. R. Douglas and N. A. Nekrasov, Rev. Mod. Phys. 73, 977 (2001) [arXiv:hep-th/0106048]. 

[19] J. Madore, Class. Quant. Grav. 9 (1992) 69. 

[20] J. GUmm and A.M. Jaffa, Phys. Rev. D 10 (1974) 536; J. Ghmm and A.M. Jaffe, and T. Spencer, 
Comm. Math. Phys. 45 (1975) 203. 



[21] W. Loinaz and R.S. Willey, Phys. Rev. D 58 (1998) 076003 [arXiv:hep-lat/9712008]; D. Lee, 
Phys. LEtt. B 439 (1998) 85 [arXiv:hep-th/9811117]. 

[22] C.-S. Chu, J. Madore and H. Steinhacker, JHEP 0108 038 (2001); B.P Dolan, D. O'Connor and 
P Presnajder, JHEP 0203 013 (2002). 

[23] B.P. Dolan, D. O'Connor, and R Presnajder, JHEP 03 (2002) 013 [arXiv:hep-th/0204219]. 

[24] N. Metropolis, A. Rosenbluth, M. Rosembluth, A. Teller, and E. Teller, J. Chem. Phys. 21, 1087 
(1953). 

[25] David Landau, Kurt Binder, "A Guide to Monte Carlo Simulations in Statistical Physics". Cam- 
bridge University Press August (2000); M. E. J. Newman, G. T. Barkema, "Monte Carlo Methods 
in statistical Physics". Oxford University Press August (2001). 

[26] J. Shao and D. Tu, "The Jackknife and Bootstrap," (Springer, New York, 1995). 

[27] William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetterling, "Numerical 
Recipes in Fortran". Cambridge University Press; 2 edition January, (1992). 

[28] Y. Shimamune, Phys. Lett. B 108 (1982) 407. 

[29] E.Brezin, C. Itzykson, G. Parisi and J. B. Zuber. Comm. Math. Phys. 59, 35 (1978); P 
Di Francesco, P. Ginsparg, J. Zinn-Justin. arXiv:hep-th/9306153; Pavel Bleher, Alexander Its. 
arXiv:math-ph/0201003. 

[30] F. G. Flores, "Simulacion de un campo escalar sobre una esfera fuzzy", PhD tesis pubUshed at 
Cinvestav, Mexico (2008). 

[31] J. Medina, W. Bietenholz and D. O'Connor, "Probing the fuzzy sphere regularisation in simula- 
tions of the 3d X^^ model," arXiv:0712.3366 [hep-th]. 

[32] Bruce M. Boghosian. Physical Review E 60 1189-1194 (1999). 



